
load("data_for_voteshares.RData")

merged <- merged2

mod2d<-glm(rhousevote~std_survey_ideology+std_survey_ideology_sq+std_leg_ideology_dems+std_leg_ideology_sq_dems+std_survey_ideology*std_leg_ideology_dems+factor(pid3)+factor(source), data=merged[merged $leg_pid3==1,],weight=weight,family=binomial(link="logit"))

mod2r<-glm(rhousevote~std_survey_ideology+std_survey_ideology_sq+std_leg_ideology_reps+std_leg_ideology_sq_reps+std_survey_ideology*std_leg_ideology_reps+factor(pid3)+factor(source),data=merged[merged $leg_pid3==3,],weight=weight,family=binomial(link="logit"))

summary(mod2d)
summary(mod2r)


# calculate the ideal point of the average Republican legislator, minus
# 0.5 standard deviations and plus 0.5 standard deviations
# "mod" stands for moderate and "ext" stands for extreme
rep_sd <- sd(merged$std_leg_ideology[merged$leg_pid3==3],na.rm=T)
rep_half_sdmod <- mean(merged$std_leg_ideology[merged$leg_pid3==3],na.rm=T) - 0.5*rep_sd
rep_half_sdext <- mean(merged$std_leg_ideology[merged$leg_pid3==3],na.rm=T) + 0.5*rep_sd
#rep_sdmod <- mean(merged$std_leg_ideology[merged$leg_pid3==3],na.rm=T) - rep_sd

# calculate the ideal point of the average Democratic legislator, plus
# 0.5 standard deviations and minus 0.5 standard deviations
dem_sd <- sd(merged$std_leg_ideology[merged$leg_pid3==1],na.rm=T)
dem_half_sdmod <- mean(merged$std_leg_ideology[merged$leg_pid3==1],na.rm=T) + 0.5*dem_sd
dem_half_sdext <- mean(merged$std_leg_ideology[merged$leg_pid3==1],na.rm=T) - 0.5*dem_sd
#dem_sdmod <- mean(merged$std_leg_ideology[merged$leg_pid3==1],na.rm=T) + dem_sd

# calculate the 5% and 95% quantiles for ideology of legislators in both parties
minrep <- quantile(merged$std_leg_ideology[merged$leg_pid3==3],.05,na.rm=T)
maxrep <- quantile(merged$std_leg_ideology[merged$leg_pid3==3],.95,na.rm=T)
mindem <- quantile(merged$std_leg_ideology[merged$leg_pid3==1],.05,na.rm=T)
maxdem <- quantile(merged$std_leg_ideology[merged$leg_pid3==1],.95,na.rm=T)

# most and least conservative districts
dists <- aggregate(merged$std_survey_ideology,by=list(merged$legislator_id),FUN="mean")
quantile(dists$x,.95)
quantile(dists$x,.05)

###
### Republicans first:
###
rdat <- merged[merged$incumb %in% c(-1,1) & merged$leg_pid3==3,]
rdat$leg_id_cd <- paste(rdat$legislator_id,rdat$leg_cd,sep="-")
gran <- 10
legislators <- unique(rdat$legislator_id)
legislators <- legislators[! is.na(legislators)]
mat <- matrix(NA,length(legislators),gran+1)
leg_pos <- rep(NA,length(legislators))
names(leg_pos) <- legislators
leg_vs <- rep(NA,length(legislators))
names(leg_vs) <- legislators
leg_vs_halfsd_mod <- rep(NA,length(legislators))
names(leg_vs_halfsd_mod) <- legislators
leg_vs_halfsd_ext <- rep(NA,length(legislators))
names(leg_vs_halfsd_ext) <- legislators
leg_vs_sd_mod <- rep(NA,length(legislators))
names(leg_vs_sd_mod) <- legislators
leg_vs_de <- rep(NA,length(legislators))
names(leg_vs_de) <- legislators
for (i in 1:length(legislators)){
    legid <- legislators[i]
    print(i)
    temp <- rdat[rdat$legislator_id==legid,]

    predicted_vs <- rep(NA,gran+1)
    pos <- seq(minrep,maxrep,(maxrep-minrep)/gran)
    # calculate predicted voteshare at legislator's actual position
    leg_vs[i] <- weighted.mean(predict(mod2r,newdata=temp,type="response"),na.rm=T,w=temp$weight)
    # legislator's actual position
    leg_pos[i] <- mean(temp$leg_ideology,na.rm=T)

    logit <- function(x){1/(1+exp(-x))}
    # calculate predicted vote shares with 1 sd below leg position, but only using direct effect
    kernel <- (predict(mod2r,newdata=temp,type="link") +
                    (-rep_sd)*coef(mod2r)["std_leg_ideology_reps"] +
                    (-2*rep_sd*temp$std_leg_ideology_reps + rep_sd^2)*
                    coef(mod2r)["std_leg_ideology_sq_reps"] +
                    (-rep_sd)*coef(mod2r)["std_survey_ideology:std_leg_ideology_reps"]
                    *.4)
                    #mean(merged$std_survey_ideology)) # the last term is effectively 0
    leg_vs_de[i] <- weighted.mean(logit(kernel),na.rm=T,w=temp$weight)               
                                  
    # calculate predicted vote shares at 1/2 sd above and below the mean party position
    temp$std_leg_ideology_reps[!is.na(temp$std_leg_ideology_reps)] <- rep_half_sdmod
    temp$std_leg_ideology_sq_reps[!is.na(temp$std_leg_ideology_sq_reps)] <- rep_half_sdmod^2
    leg_vs_halfsd_mod[i] <- weighted.mean(predict(mod2r,newdata=temp,type="response"),na.rm=T,w=temp$weight)
    temp$std_leg_ideology_reps[!is.na(temp$std_leg_ideology_reps)] <- rep_half_sdext
    temp$std_leg_ideology_sq_reps[!is.na(temp$std_leg_ideology_sq_reps)] <- rep_half_sdext^2
    leg_vs_halfsd_ext[i] <- weighted.mean(predict(mod2r,newdata=temp,type="response"),na.rm=T,w=temp$weight)
    # 1 sd below the candidate position
    temp <- rdat[rdat$legislator_id==legid,]
    temp$std_leg_ideology_reps[!is.na(temp$std_leg_ideology_reps)] <- temp$std_leg_ideology_reps[!is.na(temp$std_leg_ideology_reps)] - rep_sd
    temp$std_leg_ideology_sq_reps[!is.na(temp$std_leg_ideology_sq_reps)] <- temp$std_leg_ideology_reps[!is.na(temp$std_leg_ideology_reps)]^2
    leg_vs_sd_mod[i] <- weighted.mean(predict(mod2r,newdata=temp,type="response"),na.rm=T,w=temp$weight)
    
    # for each quantile, calculate the predicted voteshare
    for (j in 1:(gran+1)){
        temp$std_leg_ideology_reps[!is.na(temp$std_leg_ideology_reps)] <- pos[j]
        temp$std_leg_ideology_sq_reps[!is.na(temp$std_leg_ideology_sq_reps)] <- pos[j]^2
        vs <- predict(mod2r,newdata=temp,type="response")
        #mat[i,j] <- mean(vs,na.rm=T)
        mat[i,j] <- weighted.mean(vs,w=temp$weight,na.rm=T)
    }
}
rownames(mat) <- legislators
colnames(mat) <- seq(minrep,maxrep,(maxrep-minrep)/gran)
save(mat,leg_pos,leg_vs,leg_vs_halfsd_mod,leg_vs_halfsd_ext,leg_vs_sd_mod,leg_vs_de,file="rep_vote_prediction.RData")

###
### Now Democrats:
###
ddat <- merged[merged$incumb %in% c(-1,1) & merged$leg_pid3==1,]
ddat$leg_id_cd <- paste(ddat$legislator_id,ddat$leg_cd,sep="-")
gran <- 10
legislators <- unique(ddat$legislator_id)
legislators <- legislators[! is.na(legislators)]
mat <- matrix(NA,length(legislators),gran+1)
leg_pos <- rep(NA,length(legislators))
names(leg_pos) <- legislators
leg_vs <- rep(NA,length(legislators))
names(leg_vs) <- legislators
leg_vs_halfsd_mod <- rep(NA,length(legislators))
names(leg_vs_halfsd_mod) <- legislators
leg_vs_halfsd_ext <- rep(NA,length(legislators))
names(leg_vs_halfsd_ext) <- legislators
leg_vs_sd_mod <- rep(NA,length(legislators))
names(leg_vs_sd_mod) <- legislators
leg_vs_de <- rep(NA,length(legislators))
names(leg_vs_de) <- legislators
for (i in 1:length(legislators)){
    legid <- legislators[i]
    print(i)
    temp <- ddat[ddat$legislator_id==legid,]

    predicted_vs <- rep(NA,gran+1)
    pos <- seq(mindem,maxdem,(maxdem-mindem)/gran)
    leg_vs[i] <- weighted.mean(predict(mod2d,newdata=temp,type="response"),na.rm=T,w=temp$weight)
    leg_pos[i] <- mean(temp$std_leg_ideology_dems,na.rm=T)

    logit <- function(x){1/(1+exp(-x))}
    # calculate predicted vote shares with 1 sd below leg position, but only using direct effect
    kernel <- (predict(mod2d,newdata=temp,type="link") +
                    (dem_sd)*coef(mod2d)["std_leg_ideology_dems"] +
                    (2*dem_sd*temp$std_leg_ideology_dems + dem_sd^2)*
                    coef(mod2d)["std_leg_ideology_sq_dems"] +
                    (dem_sd)*coef(mod2d)["std_survey_ideology:std_leg_ideology_dems"]*
                    -.4)
                    #mean(merged$std_survey_ideology)) # the last term is effectively 0
    leg_vs_de[i] <- weighted.mean(logit(kernel),na.rm=T,w=temp$weight)   

    # calculate predicted vote shares at 1/2 sd above and below the mean party position
    temp$std_leg_ideology_dems[!is.na(temp$std_leg_ideology_dems)] <- dem_half_sdmod
    temp$std_leg_ideology_sq_dems[!is.na(temp$std_leg_ideology_sq_dems)] <- dem_half_sdmod^2
    leg_vs_halfsd_mod[i] <- weighted.mean(predict(mod2d,newdata=temp,type="response"),na.rm=T,w=temp$weight)
    temp$std_leg_ideology_dems[!is.na(temp$std_leg_ideology_dems)] <- dem_half_sdext
    temp$std_leg_ideology_sq_dems[!is.na(temp$std_leg_ideology_sq_dems)] <- dem_half_sdext^2
    leg_vs_halfsd_ext[i] <- weighted.mean(predict(mod2d,newdata=temp,type="response"),na.rm=T,w=temp$weight)
    # 1 sd above the candidate position
    temp <- ddat[ddat$legislator_id==legid,]
    temp$std_leg_ideology_dems[!is.na(temp$std_leg_ideology_dems)] <- temp$std_leg_ideology_dems[!is.na(temp$std_leg_ideology_dems)] + dem_sd
    temp$std_leg_ideology_sq_dems[!is.na(temp$std_leg_ideology_sq_dems)] <- temp$std_leg_ideology_dems[!is.na(temp$std_leg_ideology_dems)]^2
    leg_vs_sd_mod[i] <- weighted.mean(predict(mod2d,newdata=temp,type="response"),na.rm=T,w=temp$weight)
    
    for (j in 1:(gran+1)){
        temp$std_leg_ideology_dems[!is.na(temp$std_leg_ideology_dems)] <- pos[j]
        temp$std_leg_ideology_sq_dems[!is.na(temp$std_leg_ideology_sq_dems)] <- pos[j]^2
        vs <- predict(mod2d,newdata=temp,type="response")
        mat[i,j] <- weighted.mean(vs,w=temp$weight,na.rm=T)
    }
}
### SINCE THIS IS *DEMOCRATIC VOTESHARE* turn it into Rep voteshare
mat <- 1 - mat
leg_vs <- 1 - leg_vs
leg_vs_halfsd_mod <- 1- leg_vs_halfsd_mod
leg_vs_halfsd_ext <- 1- leg_vs_halfsd_ext
leg_vs_sd_mod <- 1- leg_vs_sd_mod
leg_vs_de <- 1- leg_vs_de
# names
rownames(mat) <- legislators
colnames(mat) <- seq(mindem,maxdem,(maxdem-mindem)/gran)
save(mat,leg_pos,leg_vs,leg_vs_halfsd_mod,leg_vs_halfsd_ext,leg_vs_sd_mod,leg_vs_de,file="dem_vote_prediction.RData")


# total effects
load("dem_vote_prediction.RData")
mean(leg_vs_halfsd_mod-leg_vs_halfsd_ext) # 0.9%
mean(leg_vs_sd_mod-leg_vs) #1.1%
mean(leg_vs_de-leg_vs) #
load("rep_vote_prediction.RData")
mean(leg_vs_halfsd_mod-leg_vs_halfsd_ext) # 1.1%
mean(leg_vs_sd_mod-leg_vs) # 1.7%
mean(leg_vs_de-leg_vs) #
